TempRanges

Code
spec_df<-data.frame(sci=c("Epinephelus morio","Lutjanus campechanus","Mycteroperca microlepis","Mycteroperca phenax","Rachycentron canadum","Rhomboplites aurorubens","Seriola dumerili"),
                    com=c("Red Grouper", "Red Snapper", "Gag Grouper","Scamp Grouper","Cobia","Vermilion Snapper","Greater Amberjack"))

scinm<-spec_df$sci
comnm<-spec_df$com

GOAL: We wanted a better way to examine species temperature ranges based off of real data, take ranges from a distribution to combat sampling bias, a way to semi “ground truth” general distributions, and a way to inform the depth of the area of analysis. The data we have available to us are the video surveys with related temps (bottom temp of survey?).

Tip

You can click on figures in this document to see them enlarged and zoom in!

FYI: For those interested this is possible by using the lightbox Quarto extension.

Code
#2011-2017
mosttemps<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/AllSpeciesSumcounts2011-2017.csv")
mosttemps<-mosttemps %>% dplyr::select(Year, Collection, LastOfTemp)

#2018
temps18<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/videopresenceabsence2015-2018.csv")
temps18<-temps18 %>% dplyr::select(Year, Collection, LastOfTemp)
temps18<-temps18[temps18$Year>2017,]

#All
temps<-rbind(mosttemps, temps18)
temps<-na.omit(temps)
Code
vid_ab<-read_xlsx("C:/Users/brittanytroast/Dropbox/Work/Random/Before Fixing Comp/Climate Fisheries Initiative/Data/Video_Counts_2011-2018.xlsx")

vid_ab<-vid_ab[,c(1,2,4,47:ncol(vid_ab))]
vidab_long<-pivot_longer(vid_ab, c(4:ncol(vid_ab)), names_to = "Species", values_to = "AB")

vid_spec_ab<-vidab_long[vidab_long$Species %in% scinm,]
# unique(vid_spec_abna$Species)
vid_spec_abna<-vid_spec_ab[vid_spec_ab$AB>0,]
Code
colnames(temps)[3]<-"Temp"
temps$Year<-NULL
vid_spec_abna$Collection<-as.integer(vid_spec_abna$Collection)
vid_abtemp<-left_join(vid_spec_abna, temps, by=c("Collection"))

vid_abtemp<-na.omit(vid_abtemp)
vid_abtemp<-merge(vid_abtemp, spec_df, by.x="Species", by.y="sci")

I first “squinted” at the temperature data per species to determine general curves. Here are the temperature distributions for the sum of catch at half degree temperature steps per species.

Code
#manual binning SUM
vid_bin<-vid_abtemp %>% mutate(tempbin=cut(Temp, breaks = seq(16, 31,0.5)), tempmid=as.numeric(str_sub(gsub("^(.*?),.*", "\\1", tempbin), 2))+.25)
vid_sum<-vid_bin %>% group_by(Species, tempmid) %>% summarise(sum=sum(AB))

vid_sum<-merge(vid_sum, spec_df, by.x="Species", by.y="sci")
ggplot(vid_sum, aes(x=tempmid, y=sum))+
  facet_wrap(~com, scales = "free_y", ncol=2)+
  # geom_bar(stat="identity", aes(fill=Species), color="black")+
  stat_summary(geom="bar", fun="sum", aes(fill=Species), color="black")+
  geom_smooth( method="loess", color="black")+
  scale_y_continuous(limits = c(0,NA))+
  labs(title = "Sum Abundance 0.5°C 2015-2017 SA",x="Temperature °C", y="Sum CPUE")+
  theme_bw()+theme(legend.position = "none",
                   strip.background = element_blank(),
                   strip.text = element_text(face="bold"))

Code
ggsave("AB_HalfDeg_Sum.png", width = 6, height = 6, units="in")

To address the data richness bias here are the 10-90% percentile temp ranges. This table includes these values as well as min and max temps the species was sampled at. This data is from the video trap survey data from 2011-2018 in the U.S. South Atlantic, South of Cape Hatteras, North Carolina. While it is not the whole area we are interested in, it does at least give us an idea of temperature ranges. If we had the northern data we might be able to get better cold temperature data, but luckily we care more about the upper temperature limit.

Code
sa_temp<-vid_abtemp %>% group_by(Species, com) %>% summarise(mintemp=round(min(Temp),2),ten_quant=round(quantile(Temp,prob = c(.10)),2),median=median(Temp), ninquant=round(quantile(Temp,prob = c(.90)),2),maxtemp=round(max(Temp),2))

# sa_temp %>% gt(rowname_col = "com", groupname_col = NULL) %>%
#   tab_header(title="Temperature Range °C") %>%
#   cols_label(Species="Scientific Name", mintemp="Lowest Temp",ten_quant="Lower Quantile", ninquant="Upper Quantile", maxtemp="Highest Temp") %>%
#   tab_stubhead(label="Common Name") %>%
#   tab_style(style=cell_text(style = "italic"), locations = cells_body(columns = Species)) %>%
#   tab_style(style = cell_text(color = "darkgrey", 
#         font = google_font("Source Sans Pro"), transform = "uppercase"), 
#         locations = cells_stubhead()) %>%
#   gt_theme_nytimes()


sa_temp<-sa_temp %>% mutate(comb=paste0("<b>",com, "</b><br><i>", Species,"</i>"), ten_quant=round(ten_quant, 2), ninquant=round(ninquant, 2))


sa_temp %>% gt(rowname_col = "comb", groupname_col = NULL) %>%
  cols_hide(columns=c(Species, com, median)) %>%
  tab_header(title="Temperature Range °C") %>%
  cols_label(mintemp="Lowest Temp",ten_quant="Lower Quantile", ninquant="Upper Quantile", maxtemp="Highest Temp") %>%
  tab_stubhead(label="Common Name") %>%
  tab_style(style=cell_text(style = "italic"), locations = cells_body(columns = comb)) %>%
  tab_style(style = cell_text(color = "darkgrey", 
        font = google_font("Source Sans Pro"), transform = "uppercase"), 
        locations = cells_stubhead()) %>%
  gt_theme_nytimes() %>%
  fmt_markdown()
Temperature Range °C
Common Name Lowest Temp Lower Quantile Upper Quantile Highest Temp

Red Grouper
Epinephelus morio

15.16

18.51

25.71

28.06

Red Snapper
Lutjanus campechanus

12.92

18.76

26.08

29.31

Gag Grouper
Mycteroperca microlepis

13.62

19.11

26

28.25

Scamp Grouper
Mycteroperca phenax

13.05

18.55

26.22

28.4

Cobia
Rachycentron canadum

13.56

20.02

25.63

28.66

Vermilion Snapper
Rhomboplites aurorubens

13.05

18.95

26.7

29.31

Greater Amberjack
Seriola dumerili

13.71

19.63

26.64

29.26

I further explored the data to see how many data points (each point represents one individual caught at a give temperature) we were excluding using the 10-90% range. Unlabeled text are max temperature, 90% percentile temperature, 10% percentile temperature and minimum temperature respectively from top to bottom. Labels within the boxplot show the total individuals caught (Total n), samples where individuals occurred (Sample n), and the number of individuals excluded by the 10-90% range (Out n). I figured we could use this to help us refine which range we would like to use (e.g. if we were excluding large numbers of individuals maybe we should expand our range or find natural breaks).

Code
#Single instances of abundance (for jitter)
vid_exp<-vid_abtemp %>% uncount(AB) %>% mutate(count=1)

#Add points outside of 10-90% to show how many are outside that range
temp_li<-list()
specs<-spec_df$com
for (i in 1:length(spec_df$com)) {
  dat<-vid_exp[vid_exp$com==specs[i],]
  lowup<-sa_temp[sa_temp$com==specs[i],]
  dat$inout<-ifelse(dat$Temp < lowup$ten_quant | dat$Temp > lowup$ninquant, "out", "in")
  temp_li[[i]]<-dat
 }
inout_df<-do.call("rbind", temp_li)

#Add Total n & Sampling n
ndf_1<-vid_abtemp %>% group_by(Species, com) %>% summarise(tot_n=sum(AB), samp_n=length(AB))
ndf_2<-inout_df[inout_df$inout=="out",] %>% group_by(Species, com) %>% summarise(out_n=sum(count))
n_df<-left_join(ndf_1, ndf_2)
n_df$lab<-paste0("Total n=", n_df$tot_n, "\nSample n=", n_df$samp_n, "\nOut n=", n_df$out_n)


#Boxplots for data
ggplot(sa_temp, aes(x=com))+
  geom_jitter(data=inout_df[inout_df$inout=="out",], aes(y=Temp),color="gray60", alpha=0.3, width = 0.4,height = 0.2)+
  geom_boxplot(stat="identity",aes(ymin=mintemp, lower=ten_quant, middle=median,upper=ninquant, ymax=maxtemp, fill=com))+
  geom_text(aes(label=mintemp, y=mintemp-0.5), size=3.5, fontface="bold")+
  geom_text(aes(label=maxtemp, y=maxtemp+0.5), size=3.5, fontface="bold")+
  geom_text(aes(label=ten_quant, y=ten_quant+0.5),color="white", fontface="bold",size=3.5)+
  geom_text(aes(label=ninquant, y=ninquant-0.5),color="white", fontface="bold", size=3.5)+
  geom_label(data=n_df, aes(label=lab, y=sa_temp$median), size=2.5, fontface="bold")+
  labs(title = "Video Survey Temp 10-90% Quantile" , x="Species", y="Temperature")+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+theme(legend.position = "none",
                   title = element_text(face="bold"))

Unlabeled text are max temperature, 90% percentile temperature, 10% percentile temperature and minimum temperature respectively from top to bottom. Labels within the boxplot show the total individuals caught (Total n), samples where individuals occurred (Sample n), and the number of individuals excluded by the 10-90% range (Out n).

Based on the above figure I explored a broader range of 5-95% and see if it was appropriate.

Code
#| fig-cap: Unlabeled text are max temperature, 90% percentile temperature, 10% percentile temperature and minimum temperature respectively from top to bottom. Labels within the boxplot show the total individuals caught (Total n), samples where individuals occurred (Sample n), and the number of individuals excluded by the 10-90% range (Out n). 

low_quant=0.05
high_quant=0.95

sa_temp_alt<-vid_abtemp %>% group_by(Species, com) %>% summarise(mintemp=round(min(Temp),2),ten_quant=round(quantile(Temp,prob = c(low_quant)),2),median=median(Temp), ninquant=round(quantile(Temp,prob = c(high_quant)),2),maxtemp=round(max(Temp),2))


#Single instances of abundance (for jitter)
vid_exp<-vid_abtemp %>% uncount(AB) %>% mutate(count=1)

#Add points outside of 10-90% to show how many are outside that range
temp_li<-list()
specs<-spec_df$com
for (i in 1:length(spec_df$com)) {
  dat<-vid_exp[vid_exp$com==specs[i],]
  lowup<-sa_temp_alt[sa_temp_alt$com==specs[i],]
  dat$inout<-ifelse(dat$Temp < lowup$ten_quant | dat$Temp > lowup$ninquant, "out", "in")
  temp_li[[i]]<-dat
 }
inout_df<-do.call("rbind", temp_li)

#Add Total n & Sampling n
ndf_1<-vid_abtemp %>% group_by(Species, com) %>% summarise(tot_n=sum(AB), samp_n=length(AB))
ndf_2<-inout_df[inout_df$inout=="out",] %>% group_by(Species, com) %>% summarise(out_n=sum(count))
n_df<-left_join(ndf_1, ndf_2)
n_df$lab<-paste0("Total n=", n_df$tot_n, "\nSample n=", n_df$samp_n, "\nOut n=", n_df$out_n)


#Boxplots for data
ggplot(sa_temp_alt, aes(x=com))+
  geom_jitter(data=inout_df[inout_df$inout=="out",], aes(y=Temp),color="gray60", alpha=0.3, width = 0.4,height = 0.2)+
  geom_boxplot(stat="identity",aes(ymin=mintemp, lower=ten_quant, middle=median,upper=ninquant, ymax=maxtemp, fill=com))+
  geom_text(aes(label=mintemp, y=mintemp-0.5), size=3.5, fontface="bold")+
  geom_text(aes(label=maxtemp, y=maxtemp+0.5), size=3.5, fontface="bold")+
  geom_text(aes(label=ten_quant, y=ten_quant+0.5),color="white", fontface="bold",size=3.5)+
  geom_text(aes(label=ninquant, y=ninquant-0.5),color="white", fontface="bold", size=3.5)+
  geom_label(data=n_df, aes(label=lab, y=sa_temp_alt$median), size=2.5, fontface="bold")+
  labs(title = paste0("Video Survey Temp ",low_quant*100,"-",high_quant*100,"% Quantile") , x="Species", y="Temperature")+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+theme(legend.position = "none",
                   title = element_text(face="bold"))

Lastly, I compared potential new ranges with what we had originally from a lit search.

Code
#Add old upper/lower temps to spec_df
spec_df$lowtemp=c(16,19,18.3,20.7,16.8,18.3,16.9)
spec_df$hightemp=c(25.7,29,27.1,27.5,32,27.2,29)

#Create df with old, 10-90%, and 5-95% quantile temps
all_temps<-data.frame(spec=spec_df$sci, com=spec_df$com,
                      low_old=spec_df$lowtemp, upp_old=spec_df$hightemp,
                      low_1090=sa_temp$ten_quant, upp_1090=sa_temp$ninquant,
                      low_595=sa_temp_alt$ten_quant, upp_595=sa_temp_alt$ninquant)

all_temps<-all_temps %>% mutate(comb=paste0("<b>",com, "</b><br><i>", spec,"</i>"))


all_temps %>% gt(rowname_col = "comb", groupname_col = NA) %>%
  cols_hide(columns=c(spec, com)) %>%
  tab_header(title="Temperature Range Comparisons") %>%
  tab_spanner(.,label = md("**Old**"), columns=ends_with("_old")) %>%
  tab_spanner(.,label = md("**10-90%**"), columns=ends_with("_1090")) %>%
  tab_spanner(.,label = md("**5-95%**"), columns=ends_with("_595")) %>%
  cols_label(starts_with("low")~"Lower", starts_with("upp")~"Upper") %>%
  tab_stubhead(label="Species") %>%
  tab_style(style=cell_text(style = "italic"), locations = cells_body(columns = comb)) %>%
  tab_style(style = cell_text(color = "darkgrey", 
        font = google_font("Source Sans Pro"), transform = "uppercase"), 
        locations = list(cells_stubhead(), cells_column_labels())) %>%
  tab_style(style = cell_text(align = "center", weight = "bold"),locations = cells_title()) %>%
  tab_style(
    style = cell_borders(
      sides = c("right"),
      weight = px(2),
      color = "#d3d3d3"),
    locations = cells_body(
      columns = c(upp_old, upp_1090))) %>%
  gt_theme_nytimes() %>%
  fmt_markdown()
Temperature Range Comparisons
Species Old 10-90% 5-95%
Lower Upper Lower Upper Lower Upper

Red Grouper
Epinephelus morio

16

25.7

18.51

25.71

16.82

26.95

Red Snapper
Lutjanus campechanus

19

29

18.76

26.08

17.58

27.18

Gag Grouper
Mycteroperca microlepis

18.3

27.1

19.11

26

17.24

26.71

Scamp Grouper
Mycteroperca phenax

20.7

27.5

18.55

26.22

17.37

27

Cobia
Rachycentron canadum

16.8

32

20.02

25.63

19.35

26.63

Vermilion Snapper
Rhomboplites aurorubens

18.3

27.2

18.95

26.7

17.57

27.41

Greater Amberjack
Seriola dumerili

16.9

29

19.63

26.64

18.59

27.2

Code
#Plot to visualize temp ranges easily
all_temps<-all_temps %>% pivot_longer(cols = -c(spec, com, comb)) %>% separate(name, c("temp_cat","sceno"))
all_temps$sceno<-factor(all_temps$sceno, levels = c("old","1090","595"))
temp_range<-all_temps %>% pivot_wider(names_from = temp_cat, values_from = value)
temp_range$com_nm<-gsub( " ", "\\\n",temp_range$com)
temp_range$com_nm<-factor(temp_range$com_nm,levels=c("Red\nGrouper","Red\nSnapper","Gag\nGrouper",'Scamp\nGrouper',"Cobia","Vermilion\nSnapper","Greater\nAmberjack"))


ggplot(temp_range, aes(x=com_nm, group=sceno))+
  geom_linerange(aes(ymin=low, ymax=upp, color=sceno, lty=sceno), position = position_dodge(width = -0.6), size=1.25, alpha=0.5)+
  geom_point(aes(y=low, color=sceno), position = position_dodge(width = -0.6), size=3, shape=21, fill="white", stroke=0.9)+
  geom_point(aes(y=upp, color=sceno), position = position_dodge(width = -0.6), size=3, shape=21, fill="white", stroke=0.9)+
  scale_color_manual(values = c("black", "red", "blue"), labels=c("Old", "10-90%", "5-95%"))+
  scale_linetype_manual(values = c("solid", "longdash","dashed"),labels=c("Old", "10-90%", "5-95%"))+
  labs(y="Temperature", x="Species", color="Distribution", lty="Distribution")+
  scale_y_continuous(breaks = seq(16,32, 2))+
  coord_flip()+
  scale_x_discrete(limits=rev)+
  theme_bw() + theme(axis.text.y = element_text(face="bold"),
                     legend.position = c(0.933,0.1),
                     legend.background = element_rect(color="black"))

Code
ggsave("CompareTemps.png", width=8, height=6, units="in")

Here I explore the general distribution based on the video trap survey in the southeast. I’m not sure this is the best data we could use but again it gives a general idea.

Code
#2011-2017
mostcoords<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/AllSpeciesSumcounts2011-2017.csv")
mostcoords<-mostcoords %>% dplyr::select(Year, Collection, Start_Latitude, Start_Longitude)

#2018
coords18<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/videopresenceabsence2015-2018.csv")
coords18<-coords18 %>% dplyr::select(Year, Collection, Start_Latitude, Start_Longitude)
coords18<-coords18[coords18$Year>2017,]

#All
coords<-rbind(mostcoords, coords18)
coords<-na.omit(coords)

This figure gives a smoothed heat map of density of catches. The data here are aggregated from 2011-2018. This figure is easier to use than the next one but can smooth over and hide some information.

Code
#Raw Data
vidab<-left_join(vid_spec_abna, coords)
vid_ab<-dplyr::select(vidab, year=Year, spec=Species, count=AB, lat=Start_Latitude, long=Start_Longitude)
vid_ab<-left_join(vid_ab, spec_df, by=c("spec"="sci"))

#Coastline
atl<-readOGR("C:/Users/brittanytroast/Dropbox/Work/Random/Before Fixing Comp/Climate Fisheries Initiative/Stats/CFI Hab Envelopes/Shapefiles/AtlCoastBuff.shx", verbose = F)
atl_sf<-st_as_sf(atl)
ext_atl<-extent(atl)
land<-maps::map("usa", fill=T, plot=F)
land_base<-st_as_sf(land)
land_base<-st_make_valid(land_base)
sf_use_s2(FALSE)
land_base_val<-st_make_valid(land_base)
land_atl<-st_crop(land_base_val, xmin =ext_atl@xmin, xmax=ext_atl@xmax, ymin=ext_atl@ymin-2, ymax=ext_atl@ymax)

This figure is similar to the previous however has hex bins with the number of individuals caught from 2011-2018. This smooths over less of the data.

Code
vid_ab$com<-gsub(" ", "\\\n", vid_ab$com)
ggplot(vid_ab, aes(x=long, y=lat))+
  facet_wrap(~com)+
  geom_sf(data=land_atl, inherit.aes = F)+
  stat_density_2d(aes(fill=..level..), geom="polygon")+
  labs(title = "Aggregrate Catch 2011-2018", fill="Level")+
  scale_fill_continuous(type = "viridis") +
  xlim(-82,-75) + ylim(25,36)+
  theme_bw()+ theme(strip.background = element_blank(),
                    strip.text = element_text(face="bold"),
                    title = element_text(face="bold"),
                    axis.title = element_blank(),
                    legend.position = c(0.85,0.15),
                    legend.title.align=0.5)

Code
ggplot(vid_ab, aes(x=long, y=lat))+
  facet_wrap(~com)+
  geom_sf(data=land_atl, inherit.aes = F)+
  geom_hex(bins = 20) +
  labs(title = "Aggregrate Catch 2011-2018", fill="Abundance")+
  scale_fill_continuous(type = "viridis") +
  xlim(-82,-75) + ylim(25,36)+
  theme_bw()+ theme(strip.background = element_blank(),
                    strip.text = element_text(face="bold"),
                    title = element_text(face="bold"),
                    axis.title = element_blank(),
                    legend.position = c(0.85,0.15),
                    legend.title.align=0.5)  

Here is the same heat map as above except with yearly information per species.

Code
dir.create("Spec Heat Maps")
vid_ab<-vid_ab %>% mutate(lat=round(lat, 2), long=round(long,2))

spec<-unique(vid_ab$spec)

for (i in 1:length(spec)) {
  sp_df<-vid_ab[vid_ab$spec==spec[i],]
  sp_df<-na.omit(sp_df)
  p<-ggplot(sp_df, aes(x=long, y=lat))+
  facet_wrap(~year)+
  geom_sf(data=land_atl, inherit.aes = F)+
  stat_density_2d(aes(fill=..level..), geom="polygon")+
  labs(title = unique(sp_df$com), fill="Level")+
  scale_fill_continuous(type = "viridis") +
  xlim(-82,-75) + ylim(25,36)+
  theme_bw()+ theme(strip.background = element_blank(),
                    strip.text = element_text(face="bold"),
                    title = element_text(face="bold"),
                    axis.title = element_blank(),
                    legend.position = c(0.85,0.15),
                    legend.title.align=0.5)
  filename<-paste0(getwd(),"/Spec Heat Maps/HeatMap_",substr(spec[i],1,3),"_",substr(sub(".*? ", "", spec[i]),1,3),".png")
  print(p)
# ggsave(plot = p, filename = filename, width = 5, height = 7, units = "in")

}

Here is the same hex bin map as above except with yearly information per species.

Code
dir.create("Spec Hex Maps")
vid_ab<-vid_ab %>% mutate(lat=round(lat, 2), long=round(long,2))

spec<-unique(vid_ab$spec)

for (i in 1:length(spec)) {
  sp_df<-vid_ab[vid_ab$spec==spec[i],]
  sp_df<-na.omit(sp_df)
  
  p<-ggplot(sp_df, aes(x=long, y=lat))+
  facet_wrap(~year)+
  geom_sf(data=land_atl, inherit.aes = F)+
  geom_hex(bins = 20) +
  labs(title = unique(sp_df$com), fill="Abundance")+
  scale_fill_continuous(type = "viridis") +
  xlim(-82,-75) + ylim(25,36)+
  theme_bw()+ theme(strip.background = element_blank(),
                    strip.text = element_text(face="bold"),
                    title = element_text(face="bold"),
                    axis.title = element_blank(),
                    legend.position = c(0.85,0.15),
                    legend.title.align=0.5)  
  filename<-paste0(getwd(),"/Spec Hex Maps/HexMap_",substr(spec[i],1,3),"_",substr(sub(".*? ", "", spec[i]),1,3),".png")
  print(p)
# ggsave(plot = p, filename = filename, width = 5, height = 7, units = "in")

}

We additionally discussed looking into the depth range of species so that we can customize our area of analysis to better represent an area that species may be able to inhabit if the bottom temperature is suitable.

I did a quick google/lit search of depth limits per species. Some were reputable sources like the fishing council and some were wiki, but I can dive deeper when we determine a game plan.

Code
depth_df<-data.frame(sci=c("Epinephelus morio","Lutjanus campechanus","Mycteroperca microlepis","Mycteroperca phenax","Rachycentron canadum","Rhomboplites aurorubens","Seriola dumerili"),
                    com=c("Red Grouper", "Red Snapper", "Gag Grouper","Scamp Grouper","Cobia","Vermilion Snapper","Greater Amberjack"),
                    depth_upp=c(5,9,18,0,0,18,18),
                    depth_low=c(330,189,76,100,76,122,73))

depth_df<-depth_df %>% mutate(comb=paste0("<b>",com, "</b><br><i>", sci,"</i>"))

depth_df %>% gt(rowname_col = "comb", groupname_col = NULL) %>%
  cols_hide(columns=c(sci, com)) %>%
  tab_header(title="Depth Range(m)") %>%
  cols_label(depth_upp="Upper Depth",depth_low="Lower Depth") %>%
  tab_style(style=cell_text(align="center"), locations = cells_body()) %>%
  gt_theme_nytimes() %>%
  fmt_markdown()
Depth Range(m)
Upper Depth Lower Depth

Red Grouper
Epinephelus morio

5

330

Red Snapper
Lutjanus campechanus

9

189

Gag Grouper
Mycteroperca microlepis

18

76

Scamp Grouper
Mycteroperca phenax

0

100

Cobia
Rachycentron canadum

0

76

Vermilion Snapper
Rhomboplites aurorubens

18

122

Greater Amberjack
Seriola dumerili

18

73

Here is the above data represented in the same way as the temperature ranges previously so we can compare among species.

Code
depth_df$com<-factor(depth_df$com, levels = unique(depth_df$com))
ggplot(depth_df, aes(x=com))+
  geom_linerange(aes(ymin=depth_upp, ymax=depth_low), size=1.25, alpha=0.5)+
  geom_point(aes(y=depth_upp, fill=com), size=3, shape=21, stroke=0.9)+
  geom_point(aes(y=depth_low, fill=com), size=3, shape=21, stroke=0.9)+
  labs(y="Depth (m)", x="Species", title="Depth range of Selected Species (lit search)")+
  coord_flip()+
  scale_x_discrete(limits=rev)+
  theme_bw() + theme(axis.text.y = element_text(face="bold"),
                     legend.position = "none",
                     legend.background = element_rect(color="black"))

I looked at the depth provided from the video trap survey but examining the only depth metric they have while sampling, it appears they do not sample even close to the full range of some of these species so I don’t think we can use that data.

Code
#2011-2017
mostdepth<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/AllSpeciesSumcounts2011-2017.csv")
mostdepth<-mostdepth %>% dplyr::select(Year, Collection, Start_Depth)

#2018
depth18<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/videopresenceabsence2015-2018.csv")
depth18<-depth18 %>% dplyr::select(Year, Collection, Start_Depth)
depth18<-depth18[depth18$Year>2017,]

#All
depth<-rbind(mostdepth, depth18)
depth<-na.omit(depth)


ggplot(depth, aes(x=Start_Depth))+
  geom_histogram(fill="#3788A6", color="black")+
  labs(title="Depth of Video Trap Sampling")+
  theme_bw()+ theme()

In order to refine the shape of our area of analysis, I tried to find shapefiles of the U.S. continental shelf. For some reason I could not find this easily. I did however find that around 140m is often where the cutoff occurs (shallow continental shelf). So instead I learned how to make my own shapefiles based on whatever depth I would like and found a really nice bathymetry tool/package that uses the ETOPO 2022 database hosted on the NOAA website.

Here is the bathymetric map of the U.S. Atlantic coast. Here I have added contour lines for 140m (general shallow continental shelf), 200m (second deepest range of our species), and 350m (deepest depth limit). I was surprised to find the area difference between these three depths is minimal.

Code
#| fig-cap: Depth is presented in a blue gradient from 0 to 6000 meters. Depth contour lines at 140m (light gray solid line), 200m (medium gray longdashed line), and 350m (dark grey dashed line).

#Get Bathy
bathy<-getNOAA.bathy(lon1 = -82,lon2 = -62, lat1 = 24.9, lat2 = 45)
bath_fort<-fortify.bathy(bathy)
bath_fort<-bath_fort[bath_fort$z<0,]
# bath_fort<-bath_fort[bath_fort$z> -2000,]
#Cut to 140m
# bath_cut<-bath_fort[bath_fort$z<0 & bath_fort$z> -140,]

# blue.col <- colorRampPalette(c("darkblue", "lightblue"))
# library(metR)

ggplot(bath_fort, aes(x=x, y=y, z=z))+
  geom_raster(aes(fill=z))+
  geom_contour(aes(color=c(as.factor(..level..)), lty=c(as.factor(..level..))), breaks = c(-140,-200,-350))+
  geom_sf(data=land_atl, inherit.aes = F, fill="gray70")+
  scale_fill_viridis_c(option = "G")+
  scale_color_manual(values = c("gray25", "gray45","gray65"), labels=c("350m", "200m","140m"))+
  scale_linetype_manual(values= c("dashed", "longdash", "solid"), labels=c("350m", "200m","140m"))+
  guides(fill=guide_colorbar(order=2), color=guide_legend(order=1,keyheight = 2, reverse=T), lty=guide_legend(order=1,keyheight = 2, reverse=T))+
  labs(title="Bathymetry", x=NULL, y=NULL, fill="Depth(m)", color="Contour", lty="Contour")+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))+
  coord_sf(xlim = c(-82,-62), ylim = c(24.9,45))+
  theme_bw()+theme(legend.position = c(0.833,0.134),
                   title = element_text(face = "bold"),
                   legend.box = "horizontal",
                   legend.spacing = unit(0,"mm"),
                   legend.box.background = element_rect(color="black", size=1))

Code
ggsave("Bathy_Un140.png", width =6 ,height = 6.5, units = "in")

Here I cut a shapefile to the 140m depth. This was the original depth I was going with before exploring the species depths but can customize once we agree on a depth.

Code
bath_cut<-bath_fort[bath_fort$z<0 & bath_fort$z> -140,]

#Rasterize
bath_sp<-SpatialPoints(coords = cbind(bath_cut$x, bath_cut$y))
r<-raster(ncols=297, nrows=297, xmn = -82,xmx = -62, ymn = 24.9, ymx = 45)
bath_r<-rasterize(bath_sp,r)
polys<-rasterToPolygons(bath_r)
poly_sf<-st_as_sf(polys)
poly_sf$id<-c(1:nrow(poly_sf))

#Remove grid points to isolate unwanted area
toremove<-c(10,12,14,16,9,11,13,39,41,43,40,42,44)
poly_sf<-poly_sf[!poly_sf$id %in% toremove,]

#Get Multiple Polygons to subset to 1 needed
poly_un<-st_union(poly_sf)
poly_cast<-st_cast(poly_un, "POLYGON")
poly_sf<-st_as_sf(poly_cast)
poly_fort<-fortify(poly_sf)
poly_fort$id<-c(1:nrow(poly_fort))
poly_fort<-poly_fort[poly_fort$id==4,]

#Buffer the final poly
poly<-poly_fort
# poly<-st_simplify(poly,dTolerance = .05)
# poly<-st_buffer(poly, dist = 0.1, nQuadSegs = 40)

ggplot(poly)+
  geom_sf(fill="#3788A6",color="#281E37", linewidth=0.75)+
  # geom_sf(data=poly_fort,  color="black", fill="black",alpha=0.3)+
  geom_sf(data=land_atl, inherit.aes = F, linewidth=0.5)+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))+
  theme_minimal()+theme(panel.border = element_rect(color="black",fill=NA))

Code
ggsave("PolyTest.png", height = 9, width = 7, units = "in")

I cut the bathymetric data to the under 140m depth so that we could get a better view of the depth within our possible new area of analysis. As well as showing the orginal bathymetric map with the proposed area of analysis at 140m.

Code
#crop & mask
bath_rast<-as.raster(bathy)
bath_crop<-crop(bath_rast, extent(poly))
bath_mask<-mask(bath_crop, poly)

# bath_fort<-fortify(bath_mask)
bath_df<-as.data.frame(bath_mask, xy=T)
bath_df<-na.omit(bath_df)

# bath_df<-bath_df[bath_df$layer> -150,]

ggplot(bath_df)+
  geom_raster(aes(x=x, y=y, fill=layer))+
  geom_sf(data=land_atl, inherit.aes = F, linewidth=0.5)+
  geom_sf(data=poly, color="black", fill=NA, linewidth=0.75)+
  scale_fill_viridis_c(option = "G")+
  labs(title="Bathymetry ~ < 140m", x=NULL, y=NULL, fill="Depth (m)")+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))+
  coord_sf(xlim = c(-82,-65), ylim = c(24.9,45))+
  theme(panel.border = element_rect(color="black",fill=NA),
        panel.background = element_rect(fill="white"),
        legend.position = c(0.83,0.17),
        title = element_text(face="bold", size = 13))

Code
ggsave("Bathy_CutShape.png", width = 5, height = 6.5, units = "in")

#All Bathy w/ shape outline
ggplot(bath_fort, aes(x=x, y=y, z=z))+
  geom_raster(aes(fill=z))+
  geom_sf(data=poly, inherit.aes = F, color="black", fill=NA, linewidth=0.75)+
  geom_sf(data=land_atl, inherit.aes = F, fill="gray70")+
  scale_fill_viridis_c(option = "G")+
  labs(title="Bathymetry < 140m", x=NULL, y=NULL, fill="Depth(m)")+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))+
  coord_sf(xlim = c(-82,-62), ylim = c(24.9,45))+
  theme_bw()+theme(legend.position = c(0.912,0.135),
                   title = element_text(face = "bold"))

Code
ggsave("Bathy_wShape.png", width =6 ,height = 6.5, units = "in")